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ABSTRACT 


Ho's Algorithm is reviewed and demonstrated with ana- 
lytic examples. A digital computer program is developed to 
implement the algorithm for single-input, single-output 
systems and used to identify linear continuous and station- 
ary systems which are driven with a unit step as the test 
input. Discrete realization of the continuous systems is 
obtained using the measured output-samples, to a step in- 


put, directly in the algorithm. 
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IT. REVIEW CF Ho's ALGORITHM 


A. INTRODUCTION 

Design of a control system which will perform accept- 
ably with respect to some criterion and satisfy the pos- 
sible constraints over its operating range depends on 
knowledge of the process dynamics. The dynamics might be 
linear, nonlinear, time variant, stationary or might be a 
function of environment. In any case the dynamics of the 
process must be formulated by a set of differential equa- 
tions (if it 1s a continuously operating system) or by a 
set of difference equations (if 1t is a discrete time sys- 
tem). An effective design and/or analysis can proceed 
after this formulation. 

These equations might be completely known or might be 
written down from the parameter values supplied by the 
manufacturer using the laws of physics. Sometimes only 
partial information or r:o information is supplied about 
the process; yet somehcw the dynamics must be formulated. 


This leads to the solution of an identification problem. 


Be BACKGROUND 
A number of algorithms can be found foes |- iEOu 


the solution of the identification problem. These range 


Hi 
Numbers iri brackets indicate the references at the 
end. 


from a partial identification, like the damping factor of a 
dominant complex pair of poles, to a complete identifica- 
tion of the systems dynamics using frequency and time do- 
Main techniques. 

The identification problem can be defined in general 
as: to find an internal description of the system (a set 
of differential or difference equations) from the given 
external description (input-output relation). 

Input-output relations may be defined in the time do- 
main (impulse response for continuous systems and pulse 
response for discrete systems) or in the frequency domain 
(transfer function). In the identification problem this 
relation is given in general as experimental data. 

Most of the algorithms are applicable to linear, sta- 
tionary systems only; and only this class of systems will 
be considered in this study. 

The following definitions and notations will be used. 

A continuous, linear, stationary system will be re- 


presented by 


Xe 


a oan BL 


y = Cx 


A is a nxn system matrix 


Bwisa nem @eistribukion meaieeax 


. qd 
x denotes at & 


C is a pxn observation matrix 
ieee)  loeaeenx)] “state vee ter 


= WCE) ais a emxkl dapmt vector 


iG 


y = y(t) 1s a pxl output vector 
and a discrete, linear, stationary system will be repre- 
sented by 

x(k+1) = Ay x(k) + By u(k) 


y(k) = Cy, x(k) 


(2) 


with dimensions same as (1) 
A continuous system is said to be controllable [2] if 
the matrix 


i I ] 
LB AB! A“B ) 
f 


—— << =e 


aa | (3)3 


has rank n> and is said to be observable if the matrix 


[oT! ac? : cat)%oT! wee - Fat)? ch] (4)4 


has rank n. Same definitions hold for a discrete system by 


B C. respectively. 


replacing A, B, C with A p’ <p 


D’ 
A system may always be partitioned into four possible 
subsystems [2,3] as shown in Fig. l. 
Part I: Controllable and observable 
Bawt, Li: Uncontrollable but observable 
Meobe Lit: Controllable but unobservable 


Part IV: Uncontrollable and unobservable. 


3 
Vertical dotted lines indicate partitioning. 


aa indicates transpose of C. 
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Fig.at 


As is seen from Fig. 1 only Part I gives the input-~ 
output relation. Although Part II seems to be contribut- 
ing to the output, if no energy is stored at this part 
(zero initial condition) any contribution will be the re- 
sult of a noise input only. 

If the system is of the order n and the subsystems 


are of order n,n 


T II’ n respectively then 


Loe: 


n=n+nse¢t+n +n (5) 


Identification (alternatively called realization) of 


the system in Fig. 1 will give a system of order nt: 


The above statements are given as the following 
theorem in [2]. "Knowledge of the impulse response identi- 
fies the completely controllable and completely observable 


part, and this part alone, of the dynamical system which 


generated it. This part is itself a dynamical system and 
has the smallest dimension among all realizations. More- 
over, this say eee identified by its impulse response 
uniquely up to algebraic equivalence", 

Therefore complete identification of a dynamical sys- 
tem necessitates complete controllability and observability 
of the system. Part I must constitute all the system and 
parts It, III and IV must be missing. 

In this study only completely controllable and observ- 


able systems will be considered. 


C. Ho's ALGORITHM 11,4] 

First the algorithm wiil be given for a discrete sys- 
tem as represented in (2). 

Let LS Nie be the present response at output 1 to a 
unit pulse at input j applied k time units ago. Let HL be 


the matrix 


H(K)4, AK) yg ooeeee A(K)y 


a 


(6) 


0000000808990 6 
oe0e00800@0000 8 @ 


{Ik} 
Se oooooocoeoeeo8 80 &@ @ @ JS 


Let HO » be the block matrix 


a: 2 =e oe eo e e e L 
H, H, . eee ee 
a (7) 
H 


it ee) 


Let Py, be the rank of HO Build the block matrices 


le 


HO , £=1,2,3, «.e. and each time find the rank Py. until 


P, equal Prize let this value of P, be n and £ be q;' then 


find the elementary transformation matrices P and Q such 


that 
il. OF. 6 SCO ee een 
O ai s 
d nes F 
3 ‘ : 
; O ° 
P(HO )Q =| ° i 1 ° (8) 
q e e O e 


10 


Save P and Q matrices. Let HOTS be the block matrix 


HOT = 
q 


i 
H., H. a ee reel 
: : ; (9) 
\ 


eee e ° ° e e 2q 


Then a state variable description of the system is 


found by letting 


A, = the nxn northwest corner of P(HOT,)Q 


ve) 
I 


7 the nxm northwest corner of P(HO,) (10) 
Cy = the pxn northwest corner of (HO )Q 
The following results can be stated 

— Poel =n is the minimum dimension for the 


q 
realization: 


2. By choosing P and Q suitably all possible realiza- 
tions may be obtained. 
3. Any two minimal realizations are isomorphic, that 


is there exists a nonsingular matrix W such that 


=T 
a ey 
Ds Dy 
Co woe (16) 
D5 oy 
= 
a en) 
D5 ny 


lige 


and W is given explicitly as 


=] 
= Ji Ak 
where V+ r=1,2 is given by 
! 2 -1 
a) ire ie Tot eee | 
a ce ! Ze ae) c= ! HG | Ch (a3) 
r ca 06 ae ; ne ng 


For continuous systems the algorithm follows the same 
pattern except in (6) each ES ge must be replaced by the 
impulse response and its successive derivatives evaluated 
at’ t@=80F that is, if oie) is defined as the impulse 


response between input j and output i then 


wi) 5 = re 
n(2) 34 = g(0) 5. (14) 
n(3)i4 5 §(0) 5. 


or equally ote) ie must be expanded in a power series for 


all t as 


a, Sa 2 7 lite ay aay 
1 


60 
cy uct — 
k= 
and each h(k);. must be replaced by ay. These necessitate 
that g(t) 5 must be a real analytic function. 
If the input-output relation is given in the trans- 
form domain Hs) <= for discrete systems and BU) ae EO an 


continuous systems) it must be representable as a power 


series in the form of 


12 


Bra = Daez, (s for cont. systems) (16) 


and each BGS) ae in (6) must be replaced by be after that 
the algorithm is developed in similar fashion. 

In any case Hye k = 1,2, ... constitutes a sequence of 
constant matrices. In a broad sense the identification 
problem may be restated as [iJ: Given a sequence of pxm 


constant matrices, Hye koe. LinGmkhe triple 


Ch (or A, B, C) of constant matrices such that 


Hy = CA ox 5 Kelemen, "ees (17) 


to show the mechanics of the algorithm the following ex- 
amples are presented. 
1. Example 1 Single-Input 4 Single-Output 
Discrete System (m = 1, p = 1) 


Let the Fh Kua wleace ‘os. OC 


Hy = 1 Ho = O He = -2 
H, = 6 H. = -14 ie = BO vGauirs 
Start building HO 9 Uae y- << s 
Ho, = |H, | Ss p, = 1 
2 “| 
Ay Hy ik 0 
HO. = — 7 po = 2 
Ho ion _ 7 


3 


1 O il O 
—_— lane 
0 iL 0 —*5 
O 1 
UHC.) = >D 
-2 ~3 
1 O 1 
= ———-) = 
P, (HO.) > an 
O -2 O 
1 O 
= ———* _ I i 
(HO,)Q, —7 3 ea 


14 


therefore the state equations are 


O 
x(k+1) = 

-2 
y(k) = [1 


this realization must 


Hy = Ch AS SB ’ k = dlp ee eee 
ell al 
Ba can be calculated by the Cayley-Hamilton technique: 
ii 
a “T 
k-1 
A = Q@4T+ QA = 
Dy O I D, 
- 20, some 
| at al | : 
|-2a, gay | 0 | 
where a and a, can be found as follows: 
Eigenvalues of A are 
an 
d, — a do = -2 - therefore 
k-1 _ k-1 _ ” 7 
(A,) = (-1) ae et ary =o on 
k-1l _ k-1 _ = = 
(A) = (-2) =e) sr ar. a 20, 





i ie 
| x(k) + u (k) 
-3 1 O- 
0 | x04) 


satisfy (17) 


iS 


Substituting for k = lA2eesaee 
H, = Ria Pb = (ee ao 
a 2 - (7) a = 6 
Hot 21k) 2 ea 2) ee 


which are exactly the same as the starting Hye era eer evewe 


A second set of P and Q matrices is 


which give the realization 


0 <2 i 
ik) eS X(k) + u(k) 
1 -3 0 
vik) = [12 0 | x(x) 


this realization also satisfies (17). 


For the isomorphism of two realizations, W is calculated as: 


-l 
els At 
>) 4 


We (uy Dal 


16 


By {i3): 





The two realizations have the following signal flow 


graphs: 





Fug. 2 


By Mason's Gain Formula both have the transfer function 


Z+3 


H(Z) = ToT) (ze) 


By a long division 
Hii a Pa RRO ES soe Gos | 4 ee 


the coefficients of z ~~ , k = 1,2, ... are again equal to 


Hee k"> 1,2, 


is 


2. Example 2 Two-Inputs 3» Three-Outputs 
Continuous System (m = 2, p = 3} 


Let the He k = 1,2, ... be 








1 0 1 | 
Hy = 1 2 “ Ho = A, F H, — 
2 1 -4 
—5 7 13 
Hy = —-16] , He = |16 32 
-34 38 74 
Then 
: 1 0 
HO, =) i 2 os P, a 
2 1 
1 O 0. 1 
1 2 -2 —4 
HO. =| 2 1 -3 —4 p., = 3 
0 1 f 1 
—2 —4 4 8 
-3 —4 8 14 


19 





-4 





=2ya 





therefore 


20 


A set of P 


Which give the realization 


x 
tA 


a 


— 
— 


il O 
0 O 
-2 0 
2 1 
-2 O 
—- O 


0 it 
il 1 
-4 ay 
1 0 
1 2 
2 1 


-5 


and Q matrices is 


ex 


ix 


Zi 


© 


© 


‘e 


this has the impulse response matrix 


3/2 - aie +s ene Naw eae 
G(t) ‘= pes 2e77e 
’ B/2 ee sf 5/265 Mises” + 5ea- 


Expanding the exponentials in power series and expressing 
the result as a matrix polynomial the following expression 


can be obtained. 


or — —— 


i 0 0 1. 1 iL 
G(t) = 1 2 | + | -2 -4/t +4 8 27. 
2 1 -3 -4 8 14 
heed = 7 13 
7 3 A 
| | -8 =Jls i ae] JES e2 ean e . 
-18 -34 38 74 


| eres, me 


Each element of the coefficient matrices corresponds to the 


a, 'S of (15) and the matrices themselves are exactly H'S- 


Ze 


A second set of P and Q matrices is 


1 0 0 0 0 0 
0 0 0 1 0 0 
i 0 il 0 0 0 
P, = 
2 - <a 0 i 
=" 0 2 2 1 0 
i 0 _ @- 3 £2 1 
> |e 4 ed 
3 oe 2 3 
= 
2 es Sy =F 
3 me 8-3 A 


which give the realization 


0 1 0 1 0 
x =| 0 0 if x + | 1 
0 -2 -3 1 1 
1 0 0 
‘“ = 0 ih it x 
i 0 1 


23 


Ye 


It can be shown as before that this realization also 
satisfies (17) 
To show the isomorphism of two realizations W is 


calculated as 


: 
aa T y., I 


Il 
© 
No 
= 
Fr 

| 
Ss 

| 
Es 


iE 0 al 0 0 0 
V5 = 0 1 O 1 -2 -1 
0 1 1 O -2 -3 
4Al 8 - 7 
To) eal 
(VoV. ) = 55 8 29 -16 
-7 -16 14 
3 1 1 
le 
VoVy = 8 as 8 
16 AS eS 


24 


—h 
iz T T 


25 


The W matrix also relates the two sets of P and Q 


matrices as 


W O WwW O 
25 Ee <2 meee 
0 I O I 
For this specific example 
1 O O O O 0 
0 1 O O 0 0 
Po =}1 1 ah 0 O O Py 
0 O 0 il O 0 
0. O 0 0 1 0 
! 
@ 0 0 6 0 1 
1 0 0 O 
O 1 0 0 
Qn = 2) 


IBA ANOTHER ALGORITHM 
The algorithm given in [5] for discrete systems 
(equally applicable to continuous systems) with scalar 


measurement (single-output) shows a similarity with Ho's 


Algorithm. 


26 


For the algorithm, a free dynamical system (no forc- 


ing function) represented as in (18) is considered 


x(k+1) Ay x (k) 


(18) 


y (k) Sp x (x) 


where 

A Vs) a Ween matrix 

C is a len vow wecotor 

x is a Wel “sitttte vector 

y is the scalar output 
and x(0) 1s known. Then HO, and HOT , matrices are const truc- 
ted as in Ho's Algorithm for single-input, single-output 
systems, from the measured output~sequence. An algebraic 
equivalent of the system is realized by Hetiting 


= Ih 
OT a ; 


ao 
II 


(19) 
[1 0 -«--- oO] 


Q 
[I 


Note that for single-input ,single-output systems of order 
n, q is equal to n, and op also is equal to n. Then this 
algorithm becomes a special case of Ho's Algorithm by let- 
ting 


P=1 
l (20) 


10 
II 


an 


Because a nonsingular matrix can always be reduced to an 
identity matrix by elementary transformations on its 


columns (or rows) only. 


27 


For complete identification of the system, complete 


controllability ,and observability was a necessary and suf- 
ficient condition for systems with forcing. Here, as By 


1s absent, controllability cannot be checked; instead a 
new condition called n-identification is imposed. 


A system is said to be n-identifiable if the matrix 
[ (0)! Ax (0)! AL2x(0) bees ta n-1,-(0) | (21) 
A ee ee | ea 
has rank n. Similarity of (3) and (21) is obvious for 
Single-input systems. 


Therefore it is Cees ble. to determine iN and oe tite 


the use is observable and n- -identifiable. These two 


~ 


| ee 


conditions are summed _ together and called’l-identifiability 


[5]. “ 


In fact another realization 1s possible by letting 


(ie 
HT 


-1 
A" = (HO HOG 
= (HO,) ~~ (HOT) 
; (22) 
c_ = [first row of HO_] 
D q 
This is also’a special case of Ho's Algorithm; simply 
‘ ‘ 
Pp = (HO.) > 
a (23) 
Cj Ii 
The following relation 
ie ered hk 
Ay = (Ap) (24) 


28 


holds for these two extreme cases. As HO and HOT are sym- 
metric matrices for single-input,single-output systems it 


can be written 


wig 


(Ap 


= [ (Hor) (H0,) | 


[ (x0,) ~~ J" [nor |* (25) 


i 


D 


(HO) (HOT) = A 


Example 1 is applicable to these cases just simply 


letting 


x(0) = By 


i, 


TIT. APPLICATION OF Ho's ALGORITHM 


A. IMPLEMENTATION OF Ho's ALGORITHM 

Ho's Algorithm perhaps is the one which requires min- 
imum computation among the existing identification schemes. 
Because of the simplicity of the computations it can be 
used in many applications. 

The algorithm is implemented with a digital computer 
program to see how it can be used in an off line fashion 
for identification. The results obtained can be extended 


for use in other applications. 


1. Digital Computer Programs” 


The program consists of one main program and four 
subroutines. The main program reads the data, builds the 
HO and HOT matrices and controls the subroutines. For dig- 
ital simulation of systems, data can be generated within 
the main program by replacing the data reading loop with a 
generating function. The algorithm does not stop when two 
Successive ranks of HO are equal, but uses all available 
data for additional checking. 

Subroutine SDWFD2 finds the successive derivatives 
evaluated at zero time for continuous systems; for dis- 
crete systems it is by-passed. FLAG controls this opera- 


Enon. 


5 
See pages 56-61 
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Four Newton Forward Differences are used to evaluate each 
derivative. Subroutine RANK finds the rank of HO each 
time it is constructed. Subroutine PHOQ finds the elemen- 
tary transformation matrices P and Q which put HO into its 
normal form. Subroutine MULT makes the matrix multiplica- 
tions required by the algorithm. 

The program is written for single input-output systems, 
but can be altered for multiple input-output systems; in 
fact only the reading data and building HO and HOT parts 
need changes in the main program. No change is required 


insbhiecssubmeu tines. 


pe APPLICATION TO CONTINUOUS SYSTEMS 

For discrete systems the algorithm is straight forward, 
the only limiting factor is the accuracy of the measure- 
ments. For continuous systems however some difficulties 
arise. First, as a test input signal an impulse is needed. 
For linear systems the impulse response is the derivative 
of the response to a step input. Therefore a step input can 
be used as the test input. Higher order inputs can also be 
used if it is necessary as long as the right derivatives 
are obtained. The program is written to accept the step 
response of the system. 

A second problem comes from the numerical evaluation 
of successive derivatives from the measured samples. 
Numerical evaluation of a derivative is an approximation 


even with exact measurements, when higher derivatives are 


Syl 


needed accuracy diminishes greatly. A small error in the 
input data deteriorates the results, as the numerical meth- 
ods magnify: these errors enormously. The sampling period 
is the biggest factor in this magnification as it appears 
in the denominator, and for each additional higher deriva- 
tive its power is increased by one. Magnification of error 
increases as sampling period decreases; but on the other 
hand for fast changing functions a small sampling period 

is needed. 

The examples presented in the following tables were 
Simulated digitally for unit step-inputs to the systems. 
Two parameters, DX and DD, in addition to the sampling 
period, seemed to have great effect on the accuracy of the 
obtained results. A suitable choice of the three gave re- 
peating an tie when the dimension of the HO matrix became 
larger than the order of the system. DX and DD are posi- 
mie small numbers; when finding the rank of HO, any 
number which has an absolute value less than DX is equated 
to ze by the program. When finding the P and Q matrices, 
which aut HO to its normal form, any number which has an 
absolute value less than DD is made equal to zero by the 
program. The Palaes of these parameters which gave a good 
result for this specific example are included in the 
tables. ; 

In example 3 a system with two real poles is consider- 
ed. The poles are not greatly separated and a gcod real- 


ization became possible because the step response is a 


a2 


fairly smooth function, and in the evaluation of the deri- 
vatives not too much error is involved, in fact the first 
four derivatives are almost exact. In Example 4 a system 
with three real poles, which are separated by a significant 
amount, 1s considered. In part A, the sampling rate is not 
so high. Identification of the pole at the origin is 
almost exact, the pole at -1000 is identified fairly well 
but the far pole is in error by a large amount. [In part B, 
the sampling rate is increased; while the pole at -1000 is 
identified more accurately and the far pole identification 
showed improvement, identification of the pole at the ori- 
gin began to deteriorate. In Example 5, a highly damped 
complex pair of poles gave a very accurate identification 
with a high sampling rate; as there is no Spread of poles 
a good choice of sampling rate results in an accurate real- 
iZatzen. In Example 6, a very lightly damped complex pole 
pair is considered. Again identification is very accurate 
because of the same reasons in Example 5. 

In the above examples simulation of systems was carried 
out digitally in double-precision. In general a high samp- 
ling rate with a high-precision measurement is necessary 
which makes the scheme impractical to use with measurement 
of analog signals. A different approach is necessary which 


is the subject of the next section. 


P 


Je 


EXAMPLE 3 


Transfer 

| Function 21010 
of The (s+5) (s+20) 
System , 


Sampling Per. | 0.005 
pont 


0.000001 





-3 A 
Numerically 0.3376x10 0.3000x10 -0.7500x10 


Evaluated 6 8 
Successive 0.1575x10 ale > Lec 0.6393x10 


ee ea) 1 27.000 m0. 2a carl 0.2830x107> 


Realized Popp hkeyp-dyeaes -0.1000x10° 


System ot 5 
Matrix O. LO00x10 -0.2500x10 


Realized 
Distribution 0.1000x10 
i Matrix Tf 


1 20 


0.2541x10- 


Realized 3 3 
Observation 0.3376x10 O0.3000x10 
Matrix 


Eigenvalues -0.4999 x 107 seg wml) .5.6. 
of System 5 
Matrix =—Oeebsoo x 10" + 7 0.0 
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EXAMPLE 4a 





em per wee oe 2 


| Traneter 12 
Pe MMTeST TOTT 10 
' of The s(s+1000) (s+1000000) 
| System 
Sampling Per.|0.001 
| DX ere 
DD 1.0/| largest element of HO ee ZO 
2 6 9 
Numerically 0.430x1L0 O0.9148x10 0.8746x10 
Evaluated 12 15 18 
Successive 0.8361x10 -0.7993x10 0.7641x10 
Derivatives | 
-0.7s08e10- 0,.698Ix10-° .+0.667Ga10-4 
, -0.5684x10~+3-9.1309x107? 0.2131x1077 
Realized | 
aie 0.1000x10* 0.3725x10°% -0.1991x10/ 
Matrix 
| 0.1065x107+7 0.1000x107 -0.3039x10* 


Realized 1 14 _16 
PISeribwe.on | 0+1000x10 0.3552x10 -0.1387x10 


Matrix T 


Realized : 6 
Observation 0 -430x10 0.9148x10 -0.8746x10 
Matrix 


9 


4. ; 
Eigenvalues i-0.2083 x 10° + 79 0.0 
of System 


Matrix 


_geoseu™io” 25 of 
7 


Oe lo7Uex 10 °~+ 3 0.0 
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EXAMPLE 4b 


Transfer 


Function | 1072 


of The s(s+1000) (s+1l000000) 
System 





: 
7 2 

















| -0.9625 0.1000x10 -0.9918x10 
Numerically 

Evaluated 0.8123x10'* 0.2927x101° -9.8084x107° 
Successive 





22) 29 35 


Derivatives -~0.3552x10 0.7400x10 


0.1704x10 










ly 8 S 


0.8633x10 ~“-0.5544x10 -0.5187x10 


0.1000x102 -0.1989x10~1+ -o.2083x10° 


0.6098x10~~? 0.1000x107 -0.2183x107 






Realized 
System 
Matrix 































Realized 
Distribution 
Matrix Tf 





A ie La) 


0.1000x10 -0.5963x10_ -0.2646x10 














Realized 
Observation 
Matrix 






7 2 


-0.9625 Oz Leo -0.9918x10 





-0.2083 x 10° + 4 0.0 


Eigenvalues 

of System 7 3 ; 

Matrix 0.9999 x 10° + j 0.0 
-0.2490 x 1077+ 7 0.0 
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EXAMPLE 5 


















Transfer 

Function 1000000 

of The (s+25000 + 325000) (s+25000 - 4725000) 
system 


0.000001 


Sampling Per. 
ee! 


0.000005 


1.0/| Largest element of HO | x 2.0 











6 11 






0.6764x10 0.8164x10 


& —“S6osexro-° ~~ =omlisogedlle 


3 _9.2533x102’  0.2595x107" 


~0.4927x10 








Numerically 
Evaluated 
Successive 
Derivatives 












0.1443x10 a5 











0.3408x10 

















0.9094x10~11 ~0.1250x10 


= -~0.5000x10° 


10 





Realized 
System 
Matrix 








Of OO x10 







Realized 
Distribution 
Matrix Tt 






1h 20 


0.1000x10 0.1058x10 > 












Realized 
Observation 
Matrix 






a 6 


0.6764x10 0.8164x10 





-0.2500 x 10° + 3 0.2500 x 107 


-0.2500 x 10° - 3 0.2500 x 10° 


Eigenvalues 
of System 
Matrix 
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EXAMPLE 6 


Transfer 
Function 
of The 
System 




















1001 
(s + 0.01 + 43500) (s + 0.01 - 3500) 


GOS LO al 
_ 


-0.2084x10 


Sampling Per. 


So) 1001x107 =O" 1986x107 


Numerically 
Evaluated 
Successive 
Derivatives 


Realized 
System 
Matrix 


pease 


0.2890x10~/ 


0.8198x10 71° 


0.1000x107 


14 


0.1128x10° 0.6234x10 


27 


BO 35e6x10-- omolaseo 


-0,2500x10° 


-0.1984x1072 


Realized 
Distribution 
Matrix 7 


I =, 


O00 O50 -0.2082x10 


Realized 6 


4 


Observation =10), 20s Sane as 0.1001x10 


Matrix 


Bigenvalues -0.9922 x lone + |] 0 .49998x 107 


of System _9 3 
Maier ix -0.9922 x 10 - j) 0.4999 x 10 
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Ci. DISCRETE REALIZATION OF CONTINUOUS SYSTEMS 

As is seen from the preceding examples, to find a set 
of differential equations for continuous systems is some- 
what a trial and error procedure,.as one must come up with 
a suitable triple of DX, DD and sampling period, if no 
apriori knowledge is available for the system. More impor- 
tantly, very accurate measurements are needed which is im- 
possible in an actual experiment with hardware. 

The main error source in the above procedure is the 
evaluation of derivatives numerically. If the measured 
data can be used directly, without evaluating derivatives, 
inclusion of this error to the algorithm can be. eliminated. 

This approach leads to a system which is operating con- 


tinuously but observed discretely as in Fig. 3 


T 







——» y (kT) 


u(t) Continuous ete) 


System 


rig, 3 


where T is the sampling period and u(t) is a step. If the 


system has the dynamics as in (26) 


x = Ax + Bu 


(26) 


y = Cx 
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Then 
-l 
x(t) ai | (stale x(t) + (Raa) Bu (s) } (27) 
As u(t) is a step (amplitude V) 


—] -l 
x(t) i, Wie x(t.) + {sT-A) —B y} (28) 


Defining 


o(t-t) Shs goddamn) 2} 


(29) 
_ a f-l s(sz-a) +B 
Mt-t,) ag { = } 
Substituting in 
| x(t) = O(tnt Jae y+ T(e-e,)v (30) 
letting » 
ve aT 
a | , (31) 
ee 
O 
x[(k+1)T] = 6(T)x(kT) + T(t)V (32) 
and 
y(kT) = Cx(kT) (33) 


@(T) and I(T) are constant matrices as the sampling 
period is constant, V can be replaced by v(kT) in (32) let- 
ting 


VAKT) >="8 ’ eget (0) ley tere ae (34) 


as the test input is a unit step. 
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Dynamics described by (32) and (33) can be represented 


as a discrete system of Fig. 4. 






Discrete y (kT) 





system 


Fig. 4 


If Fig. 3 is driven by a unit step and Fig. 4 with a 
unit impulse train of period T, they both will give the 
same y(kT). Therefore samples taken as the result of ex- 
periment on Fig. 3 can be used to identify the system in 
Fig. 4 as if it is driven by a unit impulse train and the 
output is measured, ‘ 

Dynamics described in (32) and (33) is an exact re- 
presentation of the continuous system of Fig. 3 at sampling 
instants and constitutes a discrete realization of the con- 
tinuous system. Two points of importance must be mentioned 
here. One is the loss of controllability and observability 
due to the sampling of continuous systems [2], as a result 
a realization which is not a faithful representation of 
the system will be obtained. If a linear-stationary—con- 
tinuous system is controllable and observable, it will re- 


main so when sampling is introduced if and only if [6] 


‘ ; k27 
Re[s, ] = Rels.,| implies Im{s;-s, J # = (35) 


4l 


s , m=1,2, ... n eigenvalues of system 


k = positives integer 


sampling period 


Ak 
If no apriori knowledge is available about the system, at 
least two realizations must be obtained with different 
sampling rates to avoid misrepresentation. 

The second point comes from the usage of a unit step 
_as the test input. As.a result the discrete equivalent in 
Piduedeits drivebpwazhtasunijt.impulise train, Ie thegouie 


ginal system has a transfer function 
H(z) = ML en ae —— (36) 


in the z-domain, the realization will be for 


Per aT ea = . teel 2 = same: 
H (z) += Geil i a,z + a5zZ + ayz + (37) 
as Yi —el,2, .6. areemeasuecd sample values to a unit 


step input. This will give an additional eigenvalue in- 
creasing the order of the system by one. This can be 


avoided by a simple manipulation of the data. 


H(z) SA) eee = asZ + (ab-a}) 2-7 + eee (38) 
meter | | 

a, = a 

an = a5 = ay (23) 

aly = a. ~ a, 
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This ‘manipulation of the data has an additional benefit, if 
there is an unknown biasing in the measurement it is remov- 
ed automatically except from the first sample. The program 
takes care of this conversion of the data. 

In the following examples systems are simulated digit- 
ally. To see the sensitivity of the algorithm to numerical 
precision, sample values were rounded-off to four signific- 
ant digits, and a realization without manipulation of the 
data is included in some examples to make comparisons pos- 
sible. As will be seen the choices of DX and DD are not 
critical as in the continuous representation of systems. 

Example 7a shows the identification of the additional 
pole at the origin because of a unit step test signal. In 
7b this pole is removed with the manipulation of the data, 
and six Significant digits are used for sample values, 
realization is very accurate. In 7c sample values are 
rounded to four significant digits, which results in a de- 
terioration in the answer. In Example 8a and 8b seven 
Significant digits are used. In 8a an additional pole is 
present because of step input, in 8b it is removed; with 
the same significant digits 8b gave a more accurate answer 
as the system order decreased by one, less computation is 
required and error is decreased. In 8c four significant 
digits are used with some deviation of the answer as a re- 


sult. In Example 9a and 9b a complex pole pair anda real 
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pole are considered, with eight and four significant digits 
of sample values respectively. Highly accurate results are 
obtained in both cases. Example 10a shows the loss of con- 
trollability and observability (as the sampling is intro- 
duced in the measurement, only observability is lost) due 


to the sampling. In 10b with a different sampling rate the 


system is identified completely. 
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EXAMPLE 7a 


Transfer 
Function 
of The 
System 


DX OgGO1 
DD 0..:0000x, 


The Type Not manipulated to take the step out 
of The Data Double-Precision 


160 
(s+1) (s+2) 





0.3552x10~" 0.0 0.7408 


1 14 


Realized 


System 0.1000x10 


1 
fiserix ~0.2464x10 


1 


0.3552x10- 


0.6661x10~1> 0.1000x101 0.2723x10 


Realized tT - 
pestribution O0.1000x10 0.3552x10 
Matias: I 


14 -0 (2220x709 © 


Realized | 
Observation 0.7244 0.2628x10 0.5374x10 
Matrix 


Eigenvalues ORS ke +o Oc 2oO x LO 
of System oat 3 
Matrix 0.9048 moO Seas— 0 


ut 


Oooo 10°«=C«~-s«C‘“}):s«éC 2777S & jigier 


Corresponding 
s-Domain 
Elgenvalues 
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EXAMPLE 7b 


Transfer 
Function 160 


of The (s+1) (s+2) 
System 


Sampling Per. Cll 
— 










Manipulated to take the step out 










The Type 


of The Data Rounded to six significant digits 






0.2220 x 1071? -0.7408 


i} 


“Realized ~ 
System 
Matrix 





















1 


QO.1000 x 10 0.1723 x ae 















Realized 
Distribution 
Matrix 


di 


A0LOOO x 10 -0.5551 x 1071® 









Realized 
Observation. 
Matrix 





at 





0.7244 0.1904 x 10 











Eigenvalues 0.8187 + 3 O80 
of System 


Matrix 













0.9048 + 43 0.0 


Corresponding -2.0 
s-Domain 
Eigenvalues -1.0 
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EXAMPLE 7c 


Transfer 
Function 
of The 
System 


= ail 


Manipulated to take the step out 















160 
(s+1) (st+2) 





The Type 
of The Data Rounded to four significant digits 


Reaitieed Olcesl x 107-” -0.7386 


System 1 ik 
Matrix 0.1000 x 10 0.17224 = LO 


Realized -16 
Dusty 1 owmt1on 0.1000 -0.6938 x 10 
Matrix T 


Realized 1 
Observation 0.7245 0.1904 x 10 
Matrix 


Eigenvalues Ore ve? 
of System 
Matrix 0.9154 


Corresponding 
s-Domain 
Eigenvalues 
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EXAMPLE 8a 


Transfer 
Function 4000 


of The ~ sas (s#1) (s+5) (s+10) 
system 


0.00001 





The Type Not manipulated to take the step out 
of The Data Rounded to seven significant digits 


7 “=—13 -13 me Te 
esis Bee 0.426x10 -0.5684x10 -0.5684x107~ -0.2030 


ees 0.1000x101 -0.5684x107> -0.5684x10"> 0.1310x10" 


0.1000x102 0. 5684x1617 -0. 2987x210" 


-0..3907x10 2" 0.0 0.1000x10! 0.2880x10! 


Reape ce : 
DarsitedsOUE DOs hal shOOUx 10 
Matrix Tf : 


us 


Rea ee aa, 
. 1 0.1077x10° 


Observation O45 27 0.2542x10 0.6169x10 
Matrix in’. 


Eigenvalues sO IS 7 0.4198 


of System 


oro aO907 3 0.2743 


O53/.L5 ) 0.2742 
0.6028 1 aQigbd 936 
Correspond- 
Ping 


s-Domain 
Eigenvalues 
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EXAMPLE 8b 


Transfer 
Function AOOO 


of The (s+l1) (s+5) (s+10) 
System 





Dales 


The Type Manipulated to take the step out 


of The Data Rounded to seven significant digits 





£5 14 


System iT 2 a 1 
Matrix 0.1000x10 0.3552x10 -0.1104x10 


0.2220x107+°? 0.1000x10* 0.1879x102 


Realized 1 
Distribution 0.1000x10 
Matrix T 


14 cS 


0.1554x10 -0.4440x10 


Realized 1 1 
Observation 0.2088x10 0.3627x10 
Matrix 


| 0.9047 + 7 0.3596 x io ** 
E1lgenvalues 
of System 


Matrix 


0.3676 - j} 0.3226 x eo 


FS 


0.6067 - j 0.3617 x 10. 


Corresponding -1.001 
, s-Domain ; 
Eigenvalues -10.003 


-4.999 
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EXAMPLE 8c 


Transfer 

Function ’ 4000 

of The . (s+1) (s+5) (s+10) 
System 





Sampling Per. Orel 


The Type Manipulated to take the step out 
of The Data Rounded to four significant digits 


14 14 


Realized | .0.1332x107 0.1776x107 0.1896 
System 1 
Matrix 0.1000x10 0.0 -0.1080x10 


0.1554x107+4  9.1000x10+ 0.1866x10~ 


1 


Realized ----}+ > -- 
Distribution 0.1000x10 
Matrix T_, 


il -14 t5 


-0.1110x10 -0.2220x10 


Realized i 1 
Observation O453 / 0.2088x10 0.3628x10 
Matrix 


14 


*Figenvalues OT 160 2Zasaalo 
of System y 
Matrix 0.1348 x #10 


14 


6.3664 mo -* 


Corresponding 
s-Domain 
Eigenvalues 
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EXAMPLE 9a 


Transfer 
Function 


of The 


80000 
(s+40) (s“+4s+400) 





of The Data 


Realized 
System 
Matrix 


Realized 
Distribution 
Matrax “L 


Realized 
Observation 
Matrix 


Eigenvalues 
of System 
Matrix 


Corresponding 
s-Domain 
Eigenvalues 


G200001 


Manipulated to take the step out 
Rounded to eight significant digits 


0.1942x1071> -0.4440x107+° 


0.1000x102 -0.4440x1071" -0,6581 


oe 0.1000x107 -0.6481 


0.1227x1072 


il -14 -14 


OS Leoox1 0 O.98332x10 0.1332580 


u 


0.4355x10 0.3682x102 1 


-0 .4645x10 
-0.3332 + j 0.7478 


-093332 - G0 7478 
-15 


oae3l1 x lo ~ + 7 ov33300Kk ae 


Mi IIB ech ae NOD 
-40.0 





el 


EXAMPLE 9b 


Transfer 
Function 
of The 
















80000 
ist40) (s+ 14544000 








fie Llype 
of The Data 


Manipulated to take the step out 
Rounded to four significant digits 









= 15 -15 -1 










-0.8881x10 -0.4440x10 OSl263x0 


| Realized 
OVS cem 
Matrix 






| 0.1000x107 0.0 -0.6576 


se 0.1000x102 














Se cp bal pt Ll gpl -0.6474 




















Realized 
Distrmoucion 
Matrix T 











Al =15 -15 


“0.1000x10 0.6661x10 0 .2220x10 












| Realized 
| Observation 
Matrix 













il al il 


On Gos Le -0.4646x10 





0.4356x10 













=0N3son = | OT 486 







Eigenvalues 
of System 
Matrix 






=03331 - 4)0n7re 


0.1885 x 107+ + 4 0.5551 x 107?° 


Corresponding 71.998 + 43 19.9 


s-Domain 


Eigenvalues pao 


=O 
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EXAMPLE 10a 
















Transfer 
Function 
of The 
System 






80000 
(st+40) (s+2+j10X) (s+2-310%) 


Sampling Per. 











The Type 
of The Data 


Manipulated to take the step out 
Double-Precision 








Realized -0.2220 
System 


Matrix 






0.1499 













0.1000 -0.8004 









Realized 
Distribution 
Matrix T 





0.1000 









Realized 
Observation 
Matrix 





OAS teglt 





Eigenvalues -0.8187 
of System 


Matrix 








Oe oor 


Corresponding -2.0 + 3 10X 
s-Domain 


| Eigenvalues 













-40.0 
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| EXAMPLE 10b 


Transfer 
Function 
of The 
System 









80000 
(s+40) (s+2+j310X) (s+2-310X) 











The Type ; Manipulated to take the step out 


of The Data =| Double-Precision 





-16 
Realized -0.2775x10 0.1108 
16 


System 0.1000x107 0.1387x107 -0.8187 
Matrix 
15 1 


0.1000x10 Og F358 


-0.8881x10+> 


--0.4440x10 


Realized . il 
Distribution | 0O.1000x10 
Matrix © 


0.1387x10 ~° 0.5551x10~ 


Realized | 1 
Observation ‘ 0.8901 0.2091x10 -0.1333 
Matrix 


0.8471 x 10°’ + 4 0.9048 


7 


Eigenvalues 
of System 


Meee 0.8471 x 10 


-~ 3 0.9048 
14 


0.1353 - j 0.9048 x 10. 


Corresponding -2.0 + 3 10% 


s-Domain 


Eigenvalues ae — 7}. LOX 


-40.0 
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ITI. CONCLUSION 


Ho's Algorithm provides a valuable tool for complete 
identification of linear—stationary systems which are com- 
pletely controllable and observable. The computations re- 
quired are simple and can be programmed easily, in fact for 
most of the computations, available library subroutines of 
a computer facility can be used. 

For discrete systems accuracy of the measurement is 
the only limiting factor. For continuous systems a dis- 
crete realization seems to be more feasible as the evalua- 
tion of derivatives introduces error. The tested examples 
showed that a five-digits=—precision data would give fairly 


good results for this approach. 
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